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Abstract 

We propose a multi-moment method for one-dimensional hyperbolic equations with smooth 
coefficient and piecewise constant coefficient. The method is entirely based on the backward 
characteristic method and uses the solution and its derivative as unknowns and cubic Her- 
mite interpolation for each computational cell. The exact update formula for solution and 
its derivative is derived and used for an efficient time integration. At points of discontinuity 
of wave speed we define a piecewise cubic Hermite interpolation based on immersed interface 
method. The method is extended to the one-dimensional Maxwell's equations with variable 
material properties. 

1 Introduction 

In this paper we develop a multi-moment method for the wave propagation, for instance, 

u t + (c(x)u) x = 0, t > 0, x e M, u(0,x) = ito(x), x£R. (1) 

Multi-moment methods approximate the solutions to differential equations using not only 
the primitive variable but also another numerical information at the grid or the cell such as 
its derivatives, the cell average of the numerical solutions. A Hermite polynomial is usually 
defined to interpolate such information. Then the numerical quantities are evolved in time 
simultaneously. 

Our method is strongly motivated by and closely related to CIP methods, which is one of the 
multi-moment methods. It was first proposed in [9] for constant velocity field. The CIP method 
uses the exact integration in time by the characteristic method and uses the cubic Hermite 
polynomial in each cell [a;,_i, xj] based on solution values and its derivatives at two endpoints 
Xj-i, xj. The method provides an accurate, less-dispersive and less-dissipative numerical so- 
lution. Here is an incomplete list of references for CIP method and related works: nonlinear 
hyperbolic equations [12], multi dimensional hyperbolic equations [11, 13], the multi-phase anal- 
ysis [15], a multi-dimensional the Maxwell's equations [7], light propagation in dielectric media 
[3], a new mesh system applicable to non-orthogonal coordinate system [14], a numerical inves- 
tigation of the stability and the accuracy [10]. The other method closely related to CIP, we refer 
to [2, 5]. 

Our contributions are as follows: Firstly, we develop the exact solution formula for solution 
and its derivative for (1) with smooth variable wave speed c(x). The cubic Hermite polynomial 
is then used to evaluate the formula approximately. Our development of the CIP scheme is 
entirely based on the characteristic method and results in a different (improved and simpler) 
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CIP scheme than the conventional one for equations with variable wave speed. For example the 
new scheme allows us to take an arbitrary time step (no CFL limitation) without losing the 
stability and accuracy. 

Next, we develop an immersed interface method [6, 16] for (1) with a piecewise constant 
wave speed. We define a piecewise cubic Hermite polynomial on a cell [xj—i, xj] which contains 
a point of jump discontinuity in c{x) using proper interface conditions and solution values and 
its derivatives at two endpoints Xj — \ ■ Xj • 

Lastly, this interface treatment is applied to the one-dimensional Maxwell's equation with 
variable material properties: We first approximate the variable coefficients by a piecewise con- 
stant (discontinuous) coefficient. The d'Alembert's based method for the Maxwell's equations is 
developed for the piecewise constant media and then applied to Maxwell system with piecewise 
constant coefficients. 

An outline of our presentation is: in Section 2 our proposed method for hyperbolic equations 
with smooth variable coefficient is developed, in Section 2.1 the error analysis of the CIP scheme 
for the constant velocity is presented, in Section 3 a CIP scheme for discontinuous media is 
developed, and in Section 4 its application to the Maxwell's system is presented. Each section 
contains some numerical tests to verify the accuracy of the proposed method. 



2 CIP method for smooth coefficient 

In this section we propose a CIP method for hyperbolic equations with smooth coefficient. We 
consider the advection equation (1) with sufficiently smooth c(x) as a model equation. The 
characteristic method is the key for designing a CIP scheme. Let u be the exact solution of 
the equation. We discretize the time domain and the spatial domain with grid size At > and 
Ax > 0. Write t n = nAt and x k = kAx, n 6 N, k G 7L. Let us consider the characteristic curve 
x = x{t) subject to x(At) = x k to the equation (1), 

= c(x(t)), x{At)=x k . (2) 

We integrate the equation backward in time to find y k = x(0). The following update formula is 
the fundamental for developing CIP scheme. 

Proposition 1. The exact solution and its derivative at time t = t n+ \, x = x k are 

u(x k , t + At) = ^^\u(y k ,t), 

c(x k ) 2 (3) 



u x {xk,t + At) 



c{x k ) c{x k ) ) c(x k ) U ^ Vk,t " > + (c(:cfc)) Ux ( yk,t ^ 



Proof. Let pi(t) = ut(t n + t,x(t)), P2(t) = u x (t n + t,x(t)) and z(t) = u(t n + t,x(t)). The system 
of ODE for pi, p2, z becomes ([4], p. 98) 

Pl {t) + c(x(t))p 2 (t) + c'{x{t))z{t) = 0, (4) 
dpi(t) 



dt 
dz(t) 
dt 



= -c'(x(t)) Pl (t), (5) 
Pl (t)+c(x{t)) P2 (t). (6) 



From (4) and (6), = —c'(x(t))z(t). Hence, from (2) 

j t (c{x{t))z(t)) = c\x(t))j t x(t)z{t) + c (x(t))j t z(t) = c(x(t)) (z(t)c'(x(t)) + j t z{i 
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Upon integrating this, we obtain z(t) = 0^ z{y k ). Similarly, from (5) pi(t) = -^^jPiiVk), 
and from (4) 

k) -MVk) + c(x(t)) P2 (x(t) + c '(x(t))4^z(y k ) = 



where pi{y k ) = -c(y k )p2(yk) ~ c'{y k )z(y k ). Hence, we obtain 

u , ( c'{y k ) c'(x(t))\ c(y k ) . ( c{y k ) \ 2 

Thus, we have (3). □ 

The formula (3) is the fundamental relations and an explicit numerical scheme is constructed 
based on the relation. Both of the primitive variable u and the spatial derivative u x will be 
simultaneously updated in CIP scheme. 

With (3) we can compute the values u(t n+ i, x k ) and u x (t n+ i,x k ) if we know the values 
u(t n ,y k ) and u x (t n ,y k ). Suppose y k falls in the j = j(k)-ih interval (xj-x, Xj). Let us denote 
the numerical solution for u(t n ,Xj) and u x (t n ,Xj) by u™ and respectively, and suppose that 
they are given at all nodes as variables. We construct the cubic Hermite polynomial on the 
interval [xj-\,Xj]; 



Hj(x) = t#_ lPl (*^i) + «X^g=i) 

+Ax t,«_ igi (^L=l) + Ax ^g 2 (^i), 



(7) 



where = (£ - 1) 2 (2£ + 1), p 2 (0 = £ 2 (3 - 2£), Ql (0 = (£ - 1) 2 £, 92(6 = £ 2 (£ - !)■ 

Approximating u(t n ,y k ) « Hj(y k ) and v(t n ,y k ) « Hj(y k ) on the interval (xj_i,Xj) and from 
(3), we arrive at CIP scheme: 



n+l _ cfafc) 
fc c(x fc ) 



n+l _ cfa) 
fe c(x fc ) 



c'(Vk) _ c'(xfc) 
c(ifc) c(x fc ) 



The cubic polynomial i/j associated with the numerical solution is called cubic interpolation 
profile on the interval [ 

For the transport equation ut + c(x)u x = 0, we have the (exact) update formula: 



In summary, CIP scheme for (1) is composed of 

CIPO Solve the characteristic ODE (2) subject to x(At) = x k to find y k = x(0). 

CIP1 For each time level t n construct the cubic interpolation profile Hj(x). 

CIP2 Update u n k +1 , v™ +1 by (8). 

Suppose that the velocity filed c is constant, there exist an integer t such that y k = x k — cAt € 
[x k -i-£, x k -i] for all k and thus the quantity A := Xk ~£ x Vk is independent of k and < A < 1, 
and the one step map (8) becomes, denoting k — £ by j, 

u n k +1 = A 2 (3 - 2A)i#_ 1 + (1 - A) 2 (l + 2A)u? 

+Ax A 2 (l - A)f™_ 1 + Ax(l - A) 2 (-A) vf, 



(10) 

+ (3A - 2)Xv^_ 1 + (1 - A)(l - 3A)^ n . 



n+l _ -6A(1-A) n 6A(1-A) n 

V k - St U j-1 A^ U j 
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We shall investigate the consistency and the stability of the scheme in the next section. 

In existing articles on CIP, the approximate scheme (8) is used in the following manner. 
One differentiates with respect to x the non-conservative form m + cu x = —c'u of the equation 
(1) to obtain m + cv x = (—c'u) x — c'u x , here we have used the notation v = u x . Next, these 
equations are split into the advection phase ut + cu x = 0, v t + cv x = 0, and the non-advection 
phase ut = —c'u, vt = (—c'u) x — c'u x . The numerical solution of the advection phase are 
given using the cubic Hermite polynomial ut = H(x*), v k = H x {x*), where the upwind point 
x* is determined appropriately for each k. Then the system of the non-advection equation 
is integrated to obtain the numerical solution of the next time level. In [12], the following 
approximate scheme for the equations is proposed: 

< +1 = (l-c'(x k )At)u* k , 

n+l _ (-, c(x k+1 )-c(x k _ 1 ) \ * K+i- u k-i- u k+i+ u k-i (11) 

J- OA™ T o A 



J V k "+ 

The other variant is proposed in [14]: 

u n k +1 = eM-c'(x k )At)u* k , v n k +1 = -u* k c"(x k )At + exp(-c' (x k )At)v* k . (12) 

Note that the CFL number is limited less than or equal to 1 in these methods. One can verify 
that (11) and (12) are approximations to (8). Indeed, if x k — y k < Ax is sufficiently small, we 
have « c %t) )At) « 1 " c '^ At « exp(-c'(x fc )At). 

2.1 Error analysis 

In this section we investigate the consistency and the stability of the CIP scheme (10) for constant 
velocity and derive an error estimates for the numerical solution. Utsumi et al [10] studies 
the stability and accuracy of the scheme through numerical simulation. We give a theoretical 
underpinning for the numerically observed evidence in terms of the eigenvalue analysis. 
The one step map (10) is written 

,n+l 



where A = — (0 < A < 1), and 

/ (1- A) 2 (1 + 2A) -(1-A) 2 A \ / A 2 (3-2A) A 2 (l - A) 

A(X) =[ , B(X) = 

\ 6A(1-A) (1-A)(1-3A) / \-6A(l-A) (3A - 2)A 

We recall that there exist an integer £ such that j = k — t for all k since the velocity is constant. 

Let us denote the exact solution of (1) for constant coefficient at the node x k at time t n by 
UJ? = u(t n ,x k ) and VJ} = u x (t n ,x k ). 

Theorem 1. 

Jjn+l \ / jjn \ / jjn \ 

A,V' ) = m ( A, V ) + B(A) ( A^, ) + °« Al > 4 >- 

Therefore, CIP scheme (10) is consistent and is of order four in space for u(x,t) and of order 
three for u x {x, t). 

Proof. This estimate directly follows from the point wise estimate for Hermite interpolation. 
Indeed from [1], denoting the Hermite polynomial of u(t n ,x) on [xj-i,Xj] by h{x), we have 



u{t n ,x)-h{x) = - ^p ,w (a;-gj_i) 2 (g-g 
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u x (t n ,x) - h x (x) = ^'^\ x - Xj_i)(aj - Xj)(x - 
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for x G [xj-\,Xj], where £o> £ij £2 £ (^j-i,^)- Thus, we have 

^ n+1 ^ = ( <*n, Vk) \ = ( KVk) + (E 1 )] 
^ Vk +1 J \ Ax Mtn,Vk) J lAx (h x (y k ) + (E 2 )™) 



Ax Vp J ' ~ v ' v V Ax VJLi J ' V Ax (E 2 )] 



where j = k — i and 



= " (4) («n,go) A 2 (1 _ A) 2 (Ax) 4 = uW(0,go-cnAt) A 2 (1 _ A) 2 

= !£^%A) A( i_ A ) (A + ^ )( Ax) 3 



n( )(0,gi-cnA f ) A(1 _ A)(A + |l )(Ax) 3. 

Here we have used the relation u^(t,x) = u^(0,x — ct). □ 

Remark 1. The consistency argument rests on the error estimates for the Hermite cubic polyno- 
mial interpolation, and the same argument can be applied to prove the consistency of the scheme 
(8) and (9); These schemes are of fourth- order accuracy in space for u and of the third-order 
accuracy in space for u x . Note that At > is chosen arbitrary independent of the mesh size 
Ax. If one takes At ~ Ax, the schemes become third-order accurate scheme in time and space 
for u and the second order in time and space for u x . 

Let us denote the error in the numerical solution by e£ = u 1 ^ — UJ} and r£ = — VZ\ From 
(10), we see that the local errors ei, r k satisfy the recursion relation 



Ax k r n k +l ) = AW { Axr] ) + BW \ Axr^ ) + { Ax{E 2 ) 



We multiply the equation by e lke and sum over k G Z; the outcome is 

e n+1 (6) \ _ ( e n {9) 



Axf n+1 {9) J °' A \Axr n (9) 



- G e,\ ( A„£nra\ )+E n (9), (13) 



where 



Ge tX = A(\)+exp(-iO)B(\), E n {9) 



El pl\n „—ik9 



^EkeziE 2 )^-' 



We investigate the stability of the map Gg )X . Let us denote the eigenvalues of the amplifica- 
tion matrix G e ,x by p li0)X and p 2 ,e,x with |pi,0,a| < \p2,e,x\- 

Theorem 2 (Conditional stability). There exists constants Mq X > andO < pg iX < 1 depending 
on (0,X) G [0, 2tt] x [0,1], such that 

\\Gl x \\ <M e , x pl x . ( 14 ) 

The proof of the theorem rests on the following observations. We postpone their proofs in 
the Appendix. 

Lemma 1. For < A < 1 and < 9 < 2ir, the eigenvalues p\ g x and p 2 9 x are simple. 
Lemma 2. \p 1AX \ < \p 2>e ,x\ < 1 for (9,X) G (0, 2vr) x (0,1). 
Now we give the proof of Theorem 2. 
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Proof. In case A = or A = 1, Gg s x becomes trivial: identity map I when A = 0, and Gg t i = 
exp(—i9)I when A = 1. Thus \\Gg J| = 1 for all n S N. If 9 = 0, 2ir, the eigenvalues of Gg t x are 
Pi,e,x = 1 — 6A + 6A 2 < 1 and P2,e,x = 1 5 an d we have 



Go,. 



1-2A 



1 



1 

1 - 6A + 6A 2 



Thus we can take pg \ = 1 and Mg \ = \\V 1 ||||^|| where V 



1-2A 



1-2A 



1 

1 



Let us assume that (9, A) G (0, 2ir) x (0, 1). The estimate (14) for Gg t \ follows from Lemma 1 
and Lemma 1. Indeed, by Lemma 1, Gg : \ is diagonalizable, i.e., there exists a invertible matrix 



Gg 



A = V a \ 



i / Pi,e,x 







(15) 



P2,0,A 

and thus ||GJ A || < \\V e ~ x l \\\\Vg. A|||p2,e,A| n for all n G N. Consequently, we have (14) by setting 
M e, P ■■= II^IIII^aII and Pe,x '■= \p2,e,x\- The inequality pg t \ < 1 follows from Lemma (1). □ 



Theorem 3. 



e n {9) 
Ax f n {9) 



< Mg,xPl, 



e°(9) 
Ax r°(9) 



'n-\ 



+ M 6 , x [Y,Pix)\\E n 



(16) 



a=0 



for all n E N. 

Proof. Apply (13) repeatedly and remember (15): 



e n (9) 
Ax f n {9) 



Gn 



e,x 



e°(0) 
Ax r°{9) 



+ (G 



n-l 



+ Gg, x + I)E n (9) 



e°(9) 
Ax f°(9) I + 



l^i=0 Pi 



0.X 



o Eto p\ a 



Vg, X E n (9), 



and so (16) follows from (14). 

It is easy to see that the residual error is of the order (Ax) 4 , i.e., 

\\E n (9)\\ ~(Ax)V 4 )(0, OU-d). 



□ 



Let us assume that the initial condition for the derivative satisfies \v^ — u x (0, Xk)\ ~ (Ax) 3 . 
Then (16) derives an error estimate: Let us denote the finial time by T and the total number of 
iteration by n. 



Corollary 1. 



\e n (9)\ ~ ^(Ax) 4 , \f»(0)\ ~ -^(Ax) 3 



for 9 G [0,2vr]. 

Remark 2. We have observed through numerical computation that 

M e,x = \\Vg~l\\\\Vg,x\\ < 3-6453, 

uniformly for (9,X) £ [0, 27r] x [0,1]. Thus, we conjecture the constant Mg t x appearing in (14) 
and (16) is uniformly bounded in 9 and A. The uniform boundness of Mg x leads to a rigorous 
proof of the error estimates 



\u(t r 



T T 

U n \\& x ~ ^- t ( Ax )\ \\u x (t n , ■) ~ V n \\ A x ~ -^( Ax f- 



We reserve the proof of the boundedness of Mg x for a future work. 
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2.2 Numerical results 



We report numerical results for the advection equation with variable coefficient. Let us consider 
the advection equation u t + (c(x)u) x = 0, c(x) = cos ^ 4 ^ +2 for x G [0, 1], t > 0, with periodic 

boundary condition. As an initial condition, we take u(0, x) = exp(— ^q q^ ); which can be 
regarded as periodic in practice. One can easily see that the solution at time t = 2, 4, 6, . . . is 
given by the initial condition u(0,x), i.e., u(2m,x) = u(0,x) for m G N. 

We report the accuracy of the CIP method. In this numerical test, the time step size is 
fixed to be At = 0.1 for each mesh size iV^ 1 , N G {50, 100, 200, 400, 800, 1600}. The numerical 
solutions at time t = 2 are computed and compared to the exact solution. The number of time 
integration is 2/ At = 20 for all N. For each mesh size, the error in the numerical solutions is 
measured by I , £ 2 and £°° norm: 



Eoo = max max \u h 
k xe[o,i] 



u(t,x k )\, 



u(t,-)\ ( 



\u(t,-)\ ei 
= {u{t,x k )Y k 



1, 2. 



(17) 



v 1 is the exact solution at 



where u n = {u k } k=1 is the numerical solution and u(t, 
the grids. The characteristic equation, dx ^ = c(x(s)), x(At) = x k , is solved backward in time 
to find the location y k = x(0) for each k. We employ Matlab build in function ode23 to solve 
the characteristic equation. 

Plot of Figure 1 (left) and Table 1 (left) show the fourth order convergence of the method. 
We also report the performance of the CIP for the transport equation ut + c(x)u x = 0, c(x) = 
cos(4ttx)+2 for x G [0,1], t > 0, with periodic boundary condition in Figure 1 (right) and Table 
1 (right). 



Advection equation 



Transport equation 




Figure 1: The fourth order convergence in space against mesh size N 1 of the numerical error 
in the numerical solution by the CIP. The error is computed by (17) at t = 2. 

Table 1: Numerical errors in the solutions for the advection and the transport equation. 





Advection equation 


Transport equation 


N 


50 


100 


200 


400 


800 


1600 


50 


100 


200 


400 


800 


1600 


ei 


1.03e-l 


1.33e-2 


6.79e-4 


5.50e-5 


2.80e-6 


2.16e-7 


1.22e-l 


1.64e-2 


8.62e-4 


7.11e-5 


3.30c-6 


2.60e-7 




9.75e-2 


1.37e-2 


7.43e-4 


6.17e-5 


3.19e-6 


2.31e-7 


1.12c-l 


1.66e-2 


9.13e-4 


7.92e-5 


3.63c-6 


2.85e-7 




1.02e-l 


1.68e-2 


9.18e-4 


9.04e-5 


4.97c-6 


4.19e-7 


l.Olc-l 


2.01e-2 


1.16e-3 


1.06e-4 


5.30c-6 


5.72e-7 



3 Immersed interface method for CIP 

In this section we develop the immersed interface method for CIP (IIM-CIP for short) for 
transport equations with discontinues coefficient; 

u t + c(x)u x = 0, t > 0, x G E, u(0,x) = u (x), x£R, (18) 
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where c > is a piecewise smooth that has a jump across x = a. We assume that the interface 
is located in an interval [xj—i,Xj]. That is, c = c~ in x < a and c = c + in x > a with 

7^ c + . We call the interval [xj—i,Xj] an irregular interval or an irregular cell, and the 
point Xj an irregular point. The interface condition on u is imposed according to the physical 
phenomena under consideration. We take [u] = as an interface condition at the interface a, 
where [u] := lim x _, a + u + (x)-lim x _, a - u~(x). Rereu~(x) = u(x)\[ Xj _ 1)(X ] andu+(x) = u(x)\ [ct)Xj] . 
The treatment of an interface condition [cu] = will be briefly discussed at the end of this section. 

From the interface condition coupled with the equation, we know the flux (cu) x = cu x is 
also continuous across the interface, and thus the derivative u x has a jump discontinuity at 
the interface. Because of the discontinuity in u x , the standard CIP using a single profile in an 
interval will not provide an accurate solution. 

We begin with the construction of a piecewise Hermite cubic polynomial that approximates 
to the solution u(t n ,x) in the interval. Next, we derive a time integration formula for the exact 
solution and its derivative. Lastly, we propose a numerical scheme (IIM-CIP scheme) to update 
the numerical solution and the derivative v™ to the next time level. 

Piecewise cubic polynomial. Let us consider to approximate the solution u(t n , x), x S [xj-x,Xj], 
The interface condition [u] = coupled with the equation (18) yields the relations [c fc |^f] = 0, 
k € N. See [16] for the derivation. Obviously, it is impossible to construct a single poly- 
nomial interpolation that satisfies the relations. We introduce two polynomials of the form 

H^{x) = Ylk=o ~w( x ~ a ) k • We determine the eight unknowns via the interface relations and 
the interpolation conditions. 

/ H~(a)=H+(a), c~ H~ (a) = c+ H+(a), 

\ {c-fH~ x (a) = ( C +)2JJ+ (a), {c'f H xxx (a) = {c+f H+ xx {a) , ^ 



H {xj-i) = uj^, 
The first equations yield 



H x {xj-i) 



H+( Xj ) 



■j ■ 



! 



c 



c a 



1 ! 



(c-fa- 2 = (c + ) 2 a+ ( C ") a a 3 - = (c + ) a a 3 + 



and thus, introducing new parameters a = (ao, • • • , a^) T , we can write H^{x) = Y^=o if 
Then using the interpolation conditions (20), we obtain the system Aa = f, where 



(20) 
(21) 

x— a \* 



.4 



/ 1 



1 



e e 2 



\ 



(e-i) 

c~ 



2(c+) 2 

(g-1) 3 
3!(c-)3 





f ul \ 




Ax v™ 










\ ^ ^-1 / 



and 9 



The determinant of A, det A 



(c-e+c+ (i-6>)) 4 



is positive for all < 6 < 1 



Ax • - Llic "^"^"""^x ^ -^i ucuxi 12(c+c-) 4 

and thus the coefficient a is uniquely determined. We define a piecewise polynomial H{x) on 
by H\[ Xj _ lt0l j = H~\[ Xj _ lt0l ], H\[ atXj ] = H + \[ aiXj ] and call it the immersed interface cubic 
polynomial to the data set /. We also denote the piecewise polynomial by H . 

Figure 2a illustrates the immersed interface cubic polynomial x) on the interval [x^—i, Xj\. 
One can see that the immersed interface cubic polynomial is continuous at the interface but has 
discontinuity in the one-sided derivative at the point. 

We have observed in Section 2.1 that the accuracy of the Hermite cubic polynomial is of the 
fourth order in the function value and is of third order in its derivative and directly affects to the 
accuracy in the one step map of the CIP method. As for the accuracy of the immersed interface 
cubic polynomial, we have the following: Let us consider the immersed interface cubic polynomial 
h (x) to the exact solution u(t n ,x). From the interface relations (19), the polynomials h + and 
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xj-1 a X j Xj-i x r c*At y, a Xj 

(a) (b) 



Figure 2: (a) The piecewise cubic polynomial on the irregular interval. It is continuous at the 
interface a and has jump discontinuity in the one-sided derivative at the interface, (b) Graphical 
illustration of the IIM-CIP scheme. The value H + (xj — c + At) can be used for updating the 
numerical solution at the time level t n+ \ at the grid Xj. 



hr can be written in the form h+(x) = £Lo ff and h~(x) = £f =0 f {^tY ■ Thc 

interpolation condition with exact solution as a data set, 

h 1) — u(t n: Xj— i), h x (xj—\) — u x (t n , -Ej— i)j h {xj) — u{t nj Xj^)j h x {xj) — Ux(t n: Xj^) : 

leads to the system Aa = f for the coefficient a, where 

f = (u + (t n ,Xj), Axu+(t n ,x), u~(t n ,Xj-i), Ax u~(t n ,x)) T . 

Theorem 4 (Accuracy of the immersed interface cubic polynomial). Let u(t,x) be the solution 
of (9) . We have the following error estimates 

u(t n ,x) - h(x) = 0((Ax) A ), u x (t n ,x) - h x (x) = 0((Ax) 3 ), x G [xj-i, xj]. 

Here the one-sided derivative is taken at the interface a. 

Proof. The Taylor expansion for u^(t n ,x) at the interface a and the interface relations 

u~(t n ,a) = u + (t n ,a), c~u~(t n ,a) = c + u^(t n ,a), 

( c ) u xx(trn o>) = (c + ) u xx (t n ,a), (c ) u xxx (t n ,(x) = (c + ) u xxx (t n ,a), 

lead to 

u+(t n , Xj ) = 6 + 6 1 ^ + |(^) 2 + |(^) 3 + O((Ax) 4 ), 
Axu+(t n ,x) = ^ + 6 2l -| F + | I ^ F + ((Ax) 4 ), 

trC^-i) = 6 + 6i^ + |(^ 1 ) 2 + |(^ 1 ) 3 + O((Ax) 4 ), 
Axu-(t n ,x) = ^ + b 2 pjz + b fpf + 0((Ax)*), 

where &o = u + (t n ,a), b\ = c + u^(t n ,a), b 2 = (c + ) 2 u^ x (t n ,a) and 63 = (c + ) v£ xx (tm a). Thus 
we have A(a — b) = (n, r2,r^,r4) T where = 0((Ax) 4 ) for i = 1,2, 3, 4. 

Since the components of A is of order O(l), the components of the inverse A -1 is also of order 
O(l). Therefore, the equation A(a — b) = r implies that a^ — b^ = 0((Ax) 4 ) for k = 0, 1,2,3. 
Thus we obtain the desired estimates. □ 
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Update formula. Now that we have constructed the immersed interface cubic polynomial in the 
irregular interval, we attempt to build an update formula for numerical solutions with the aid 
of the method of characteristic. Throughout of this section we assume that CFL number is less 
or equal to 1, i.e., c+At < Ax. Thus yk is always included in the interval [xk-i,Xk] for all k. 

If Xj is an regular point, the numerical solution is updated by (10) with A = c ^ t or A = C A ^ 
depending on the velocity at the grid x = Xj. Let us assume that Xj is an irregular grid, and let 
us consider the characteristic curve x = x(s) with x(At) = Xj. Let us denote the upwind point 
z(0) by yj. 

Proposition 2. Suppose that there exists an interface in the interval [xj-±, Xj]. If a < Xj— c + At, 
Uj = Xj — c+At, 

u + (t n+ i,Xj) = u + (t n ,yj), u+(t n+1 ,Xj) = u+(t n ,yj). 

Ifxj - c+At < a, 

yj = a+ ^f(xj - a) - c~At, 

u+(t n+1 ,Xj) = u-(t n ,yj), u+{t n+ i,Xj) = ^pu~(tn,yj). 

Proof. Let us assume that a < Xj — c + At. The characteristic ODEs (2)— (5) become x(s) = c+, 
z(s) = p'i(s) = p2(s) = 0, and thus we have yj = Xj — c + At, u + (t n+ \,Xj) = u + (t n ,Xj — c + At) 
and u + (t n+ \,Xj) = u+(t n ,Xj — c + At). 

Next let us assume that Xj — c + At < a. There exists s* > such that a = x(s*). The 
characteristic curve obeys 

x(s) = c~, (0<s<s*), i(s) = c + , (s* < s < At). 

It is obvious to see that a — yj = c~s*, Xj — a = c + (At — s*), which yields yj = a + 
^f(xj — a) — At, s* = At — X] c+ a ' . The characteristic ODE for z{s) = u(s,x(s)) becomes 
£u~(s,x(s)) = 0, (0 < s < s*) and -^u + (s,x(s)) = 0, (s* < s < At). Using the interface 
condition [u] = 0, it follows u + (t n+ i, Xj) = u + (s*,a) = u~(s*,a) = u~(t n ,yj). Similarly, with 
the characteristic ODE for u x (s,x(s)) and the relation [cu x ] = 0, we obtain u x (t n+ i,Xj) = 
u+(s*,a) = ^u~(s*,a) = ^u~{t n ,yj). □ 

Thus we arrive at a CIP scheme for the equation (18): Let H + and if" be cubic polynomials 
defined by (19) and (20). If a < x d - c+At, 



u 



JJ+ 1 = H+{ X j - c+At), v] +l = H+{ X j - c+At). (22) 



and, if Xj — c+At < a, 



u 



^ = H~( yj ), v ^ = ^- H -{yj). (23) 



where yj = a + fr( x i — ct) — c^ At. 

It seems one must appropriately choose either (22) or (23) according to the condition a < 
Xj — c+At or Xj — c+At < a. However, no need for this special treatment arises and the update 
scheme at the irregular point is solely based on (22). Indeed we can show that if Xj — c+At < a, 
the numerical solution is also given by (22), i.e., 

u n+1 = H-(yj) = H+( Xj - c+At), v n+l = ^ff-foi) = H+( Xj - c+At). (24) 
From ^r(yj — a) = Xj — c+At — a and (21), we have 

/c+V 

a e (yj - a) 1 = a\ ( — J (yj - a) 1 = aj(xj - c+At - a) e , 
c~ /c+A^ 1 

— aj(yj - a) 1 ' 1 = a\ ( — j (yj - a) 1 ' 1 = aj(xj - c+At - af" 1 , 
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which implies (24). 

Figure 2b illustrates that we can use H + (xj — c + At) and its derivative to update the numerical 
solutions and v r - +1 at the irregular grid even when Xj — c + At < a and (23) is not necessary 
for the computation of the numerical solution. 
Here is the summary of the IIM-CIP. 

IIM-CIP1 Construct the cubic profile on each interval. If the interval is irregular, then con- 
struct the immersed interface cubic polynomial. 

IIM-CIP2 Update u™ +1 and If Xj is a regular grid, use (10) with A = ^ where c = c~ 

if Xj < a or c = c + if a < Xj. If Xj is the irregular grid, use (22). 

It is worth pointing out that if c is a continuous constant, then H + (x) = H~(x) by (21), 
and thus the immersed interface cubic polynomial H(x) is identical to the cubic interpolated 
polynomial (7). 

Let us consider the advection equation Ut — {cu) x — 0, with c > 0. Assume that an interface 
locates x = a G [xj-i,xj]. Then Xj-i is an irregular point and the numerical solutions at Xj-i 
are give by 

u]+l = H~ { Xj - 1 + c~ At ) , v^+l = H~ (xj^ + c-At). 
Interface condition [cu] = 0. 

The interface condition yields to the interface relations [c l ^f] = 0, i £ N. We define the 
immersed interface polynomial using the conditions 

[cu] = 0, [c 2 u x ] = 0, [c 3 u xx ] = 0, [c 4 u xxx ] = 0, 

and the interpolation condition (20). The numerical solutions are given by the same form as 

u n+l = H +( x . _ C+At ^ y n+l = H +( x . _ C+At) _ 

We list up the features and the advantages of IIM-CIP: 

(i) The method becomes the standard CIP if the discontinuities in the coefficients disappear. 

(ii) The structure of the IIM-CIP is as simple as CIP; The cubic interpolation profile H is 

replaced by the immersed interface cubic polynomial. 

(iii) The immersed interface cubic polynomial is an interpolation for the solution on the irregular 
cell, and the order of accuracy is four for u(t, x) and three for its derivative. Therefore, 
the order of accuracy is maintained to be the same as that of the standard CIP scheme 
for CFL < 1. 

(iv) No grid refinement is required to maintain the accuracy in a vicinity of the interface. 



3.1 Numerical results 

In this section we present some numerical results to illustrate IIM-CIP for discontinuous velocity 
and verify third order accuracy in time and space. 

Example 3.1.1. We consider the transport equation (18) with discontinuity 



c(x) 



c\ < x < a, 
C2 a < x < 1. 



The periodic boundary condition is imposed. The interface condition [u] = and [cu] = 
at the interface x = a are tested. We apply IIM-CIP method developed in section 3 to this 
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problem. The constants a, c\ and C2 are given below. We take u(0, x) = exp(— ^ q 5 ^ ), as 
initial condition. The exact solution for [u] = is given by u(t,x) = u(0,y), where 

{x — C2t, if x > a and t < 

x — c\t, if x < a, 

^(x — C2t) + (1 — otherwise . 

The exact solution for [cu] = is u(t,x) = ^||yit(0,y). 

In the numerical experiments, the spatial domain [0, 1] is uniformly discretized with mesh 
size Ax = i for N G {50, 100,200,400,800, 1600}, and the time step size is At = 0.5 Ax. The 
location of the interface is set to be a = 0.5. The velocity is set to be c\ = 1 and C2 = 2. 
The initial condition for the numerical solution {u^}^ is given by u Q k = u(0,Xk), and {v^}^ is 
computed by central finite difference of {u®} 1 ^. For each mesh size, the error in the numerical 
solutions at time t = 0.4 is measured by (17). 

Plots of Figure 4 show numerical solutions (dot) and exact solutions (solid line) at time 
t = 0.14 (left), t = 0.17 (center) and t = 0.25 (right) to the transport equation with jump 
condition [u] = at a = 0.5. The mesh size ^gg is used to compute the numerical solutions. The 
vertical solid line indicates the location of the interface. As the wave passes the interface, it slows 
down and becomes narrower. No spurious oscillation is observed in the numerical solutions. 

Numerical solutions and the exact solution for the interface condition [cu] is plotted in Figure 
5. The exact solution u has a jump discontinuity at the interface. We see that the numerical 
solution also exhibits the distinct jump at the interface. A magnification of the figure at time 
t = 0.14 around the interface is depicted in Figure 6 so that the jump discontinuity in the 
numerical solution and the exact solution are more visible. We observe that the numerical 
solution almost coincides with the exact solution. 

Plots of Figure 3 show the error in the computed solutions at t = 0.4 against mesh size iV -1 . 
Grid refinement studies confirm that the third order accuracy in time and space is achieved at 
all grid points. The third-order accuracy in time and space of the standard CIP is maintained 
even in the presence of the interface. 

In [16], an immersed interface method is presented for a piecewise constant velocity. A 
piecewise quadratic interpolation on a cell i, Xj] which contains a jump discontinuity in 
c(x) is constructed based on the immersed interface method and solution values at local three 
grid points. The underlined time integration method used in [16] is Lax-Wendorff method. The 
method causes visible oscillations in the numerical solution due to the fact that the Lax-Wendroff 
scheme is dispersive. 



[u]=0 [cu]=0 

iu! , ., 10° t , — — 




Figure 3: The third order convergence against mesh size N of the numerical error in the 
numerical solution by the IIM-CIP. The error is computed by (17) at t = 0.4. 
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Figure 4: 1-D transport with c~ = 2 and c + = 1 and jump condition [u] = 0. The vertical solid 
line indicates the interface a = 0.5. The plots are the numerical solution (dot) and the exact 
solution (solid) at t = 0.14 (left), t = 0.17 (center), t = 0.25 (right). The mesh size is ^gg. 




Figure 5: 1-D transport with c~ = 2 and c + = 1 and jump condition [cu] = 0. The vertical solid 
line indicates the interface a = 0.5. The plots are the numerical solution (dot) and the exact 
solution (solid) at t = 0.14 (left), t = 0.17 (center), t = 0.25 (right). The mesh size is ^ . 



4 IIM-CIP for Maxwell's equations in one dimension 

Let us consider one dimensional Maxwell's equations 

(25) 



eE t — H x 



uH t — E x 

for x G [0,1], t > 0, with periodic boundary condition. 

We begin with the case u and e are constants. The solution update scheme is based on the 
D'Alambert formula: 

tt(. \ _ H{t n ,x-cAt)+H(t n ,x+cAt) _ E(t n ,x-cAt)-E(tn,x+cAt) 

n \ l n+l, X) — 2 2c fi J 

p /. \ _ E{t n ,x-cAt)+E(t n ,x+cAt) _ H{t n ,x-cAt)-H(t n ,x+cAt) 

I - J \J"n+li X) — 2 2ce ' 

where c = The formula follows from the fact that the transformations 



- ^~iH - yfeE, u 2 = ^JIH + y/eE, 



u 



satisfy the transport equations 



u\ + cu x = 0, uf + cu 2 = 0. 



Let us denote the numerical solution to H(t n ,x k ), H x (t n ,x k ) E(t n ,x k ) and E x (t n ,x k ) by 
H%, DH£, E% and DE^, respectively. Let h k ~ 1)k {x) and e k ~ l,k (x) be the cubic polynomials on 
the grid [x k -i,x k ] at time t = t n determined by the conditions 

/ i fc " 1 ' fc (x A: _ 1 ) = H%_ v ^- 1 '*(x fc _i) = DEl_ x h k -^ k (x k ) = HI ht~ hk (x k ) = DH%, 
e k - x ' k {x k ^) = E%_ v e*" 1 '*^-!) = DEl_ x e k ^ k {x k ) = E% e^ 1,fe (x fe ) = DE£. 

Let us assume the CFL number is less or equal to 1. Then x k — cAt £ [x k —i, x k ] and CIP scheme 
for equation (25) is given by 

Tjn+l _ h k - 1 ' k (x k ~cAt)+h k ' k+1 (x k +cAt) _ e k - 1 ' k (x k -cAt)-e k ' k+1 (x k +cAt) 
n k ~ 2 ~~ 2cfi ' 

pn+l _ e k ~ 1 - k (x k -cAt)+e k ' k+1 {x k +cAt) _ h k - 1 ' k (x k -cAt)-h k ' k+1 (x k +cAt) 
~ 2 2ce ' 
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Figure 6: A magnification of the left plot in figure 5. The plots are the numerical solution (dot) 
and the exact solution (solid) at t = 0.14. 



and DH^ +1 and DE'^ +1 are given by replacing the polynomials within H r £ +l and E^ +1 with 
the derivatives. 

Now we consider the Maxwell's equations with e and [i being piecewise constants, i.e., e = e~ 
and \x = fj,~ in x < a and e = e + and ^ = /i + in a < x. The interface conditions for H and E 
at the interface are [H] = and [E] =0. The interface relations are 



[H] 
[E] 













1 


o, 


-H x 


= 0, 


— H X x 


= 0, 






e 




lie 
















1 


o, 


—E x 


= 0, 


^xx 


= o, 






_V _ 




lie 




jx 2 e 



-E^r 



0, 



which are obtained in a usual manner in IIM; first differentiating the relation [H] = with 
respect to t, then substituting the equation fxHt = E x to get jE x 



0. Differentiate this 



0. The others 



relation again with respect to t and substitute eEf = H x , we obtain — H x 
are given by repeating this procedure. 

Let us assume that the interface x = a is included in [xj-\,Xj]. Let us denote the immersed 
interface cubic polynomials for H and E on the cell [xj_i,Xj] by h^(x) and e ± (x) respectively. 
From the interface relations, 

h (x)-a + e ai — + {ne) y(^) + (/*) gfC"^-) » 



-l- 63 x — a, 



e±(x) = b + + {,e)^C-^? + (A)^(^) 3 . 

Ax 2 Ax 3! Ax 

The coefficients a and 6 are determined by the interpolation condition at two end points Xj, 
i.e., 

A(e, /i)a = {Hf, Ax DiZ?, H^_ v Ax DH*_ x ) r , 



A(n, e)b = (£™, Ax DEJ, E%_ lt Ax DE™_ X 



where 



A(e,fi) 



(l e+9 ^f* \ 

e + (ne)+e b*p> 

1 e - (fl _i) (HtM (gdqgzig 

2 / 



Xj — a 
Ax 



\0 e~ (fie)- (9-1) 
The numerical solutions H™ +1 and at the irregular point x 7 - are then given by 



H n + 1 



n+1 



h+ (xj -c+ At)+hi-i +1 (XJ+C+ At) _ e+(x j -c+At)-ei>i+ 1 (x j +c+At) 
2 ~~ 2c+fi+ ■ 

e+ (xj-c+ At)+ei'i +1 (xj+c+ At) _ h+ (xj -c+ At)-h^i +1 (xj +c+ At) 
2 2c+e+ ' 



(26) 
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where c + = , 1 DH™ +1 and DE™ +1 are given by replacing the polynomials within H™ +1 

and E" +1 with the derivatives. We note that the point Xj—i is also an irregular point, and the 
immersed interface cubic polynomials are used to update the numerical solutions at this point. 
For instance, H™_i is given as 

TTri+l _ h?- 2 'i- 1 (xj-i—c~ At)+h" (xj-!+c~ At) e J_2j ~ 1 (x : ,_i-c~Ai)-e-(xj_i+c~Ai) 
j-l ~ 2 2c- /x- ' 

where c" 



(i e 



4.1 Application to Maxwell's equations with variable material parameters 

We apply the method we developed for variable e(x), y(x). We approximate e(x), /i(x) by the 
piecewise constant (discontinuous) media: 

1 f x j^"~T _ , r x i + ^r 

£ ( x ) = Ai / e ( x ) dx ' ^\ X ) = 35" / M(^) dx, 

/ ry> . Ax It-. Ax 

J X 3 — J X 3 — 

on (xj-x/2, x j+i/2) f° r au 3- Then one can apply (26) for the both backward and forward manner 
to obtain 

o-n+l _ rf^jixj-cA^+h-j^Xj+cAt) _ e+^^Xj-cA^-eJ^^Xj+cAt) 

n j - 2 _ 2c/j,_ ' (27) 

n+1 _ et^.^-cAQ+e^.^Orj+cAt) _ ^ ^ (xj - c At) - j+ 1 ( Xj +c At ) ^ > 

^3 ~ 2 2ce • 

By replacing the polynomials within (27) with the derivatives, we obtain DHJ +1 and DE? +1 . 
Here c = 5j, y = ftj, e = Ej on the cell (^-1/2)^+1/2)5 an d j( x )> e f-i j( x ) are the 
immersed interface cubic polynomials on [xj-i,Xj] and hj- +1 (x), e- +1 (a;) are the immersed 
interface cubic polynomials on [xj, £Ej+i]. 

4.2 Numerical results 

We present two examples to illustrate the potential of the IIM-CIP for the Maxwell's equations. 
Example 4.2.1. Consider the Maxwell's equations (25) with 

_ f e~ = 1, < x < a, _ f yr = 1, < x < a, 

£W - \ e + = a < x < 1. ' MXj ~ \ M + = 3, a < x < 1. 

The spatial domain [0, 1] is uniformly discretized with mesh size Ax = ^gg, and the time step 
size is At = 0.5Ax. The location of the interface is set to be a = 0.5. As an initial condition, 

we take H(0,x) = exp(- {x ~^f ) and E(0,x) = -J^H(0,x). 

Plots of Figure 8 show the numerical solutions to H (t, x) (left column) and E(t, x) (right 
column) at time t = 0, t = 0.3, t = 0.35 and t = 0.5. There are no spurious oscillations observed 
in the vicinity of the interface, at least for this example. 

Example 4.2.2. Consider the Maxwell's equations e(x)Et = H x , n(x)Ht = E x for x G [0, 1], t > 
with periodic boundary condition, where e(x) = /i(x) = \ cos(4-7rx) + l. As an initial condition, 

we take H(0,x) = exp(— ^q 5 ^ ) and E(0,x) = 0. We apply the IIM-CIP developed in Section 

4.1. In this numerical test, the time step size is chosen to be At = — — = 0.25iV -1 for 

each N G {50, 100, 200, 400, 800, 1600}. The numerical solution is integrated in time by (27). 

The numerical solutions of the magnetic field H at time t = 1 are compared to the exact 
solution, which is identical to the initial condition, i.e., H(l,x) = H(0,x). For each mesh size, 
the error in the numerical solutions is measured by i 1 , £ 2 and £°° norm: 

eoo = max max \H% - H(l, x k )\, * = 1 , - j, \, , i = 1, 2. (28) 
k x£[0,l] 
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Figure 7 shows errors in the numerical solutions against mesh size N . Grid refinement studies 
confirm that the second-order convergence in time and space is achieved. The second-order 
accuracy in the approximation of e and p, results in the second-order convergence in the numerical 
solutions. 




Figure 7: The error (28) in the numerical solutions at t = 1 of the Maxwell's equation with 
variable material parameters e(x) = /jl(x) = \ cos(47rx) + 1 computed by IIM-CIP against mesh 
size iV -1 . The second order convergence in time and space is observed. 



5 Conclusion 

We have developed a numerical scheme for one-dimensional hyperbolic equations with variable 
coefficient. The method is based on the backward characteristic method and uses the solution 
and its derivative as unknowns and cubic Hermite interpolation for each computational cell. 
The consistency and the conditional stability of the method was presented. We have proposed a 
numerical scheme for one-dimensional hyperbolic equations in a discontinuous media. We have 
constructed the immersed interface cubic polynomial. We have extended the method to the 
one-dimensional Maxwell's equations with variable material properties by approximating with 
a piecewise constant media. 
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7 Appendix 

Proof of Lemma 1 

Let p{z) be the characteristic polynomial of Gq^x. p(z) is written as 

p(z) = z 2 + 2 (e- ie \ (1 - 3A + A 2 ) - (l - 2A + A 3 )) z + (1 - A) 4 + e" 2ie A 4 - 2e" 4e A (l - 2A 2 + A 3 ) 
= z 2 + (3z + 7. 

The necessary and sufficient condition for pi t e,x 7^ P2,e,x is that (3 2 — 47 7^ 0. We obtain, after 
some work, 

/3 2 -4 7 = 4(-l + A) 2 A 2 (cos#-isin#)[2 (5 + A - A 2 ) + (-1 - 2A + 2A 2 ) cos0 + 3i(-l + 2A)sin0], 
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and so j3 2 — 4 7 / is equivalent to 

q(X, 9) := (2 (5 + A - A 2 ) + (-1 - 2A + 2A 2 ) cos Of + (3(-l + 2A) sin 9) 2 / 0. 

The last term (-1 + 2A) sin (9 equals to when A = 1/2 or sin 9 = 0. But g(l/2,0) = 3(7 - 
cos 6>)/2 / 0, g(A,0) = 9 and g(A,7r) = 11 + 4A - 4A 2 . Thus we see that q(X,9) ^ for all 

< A < 1. □ 
Proof of Lemma 2. 

We employ the theory of Schur to check whether the roots of the polynomial p(z) reside inside 
the unit circle. Let p*(z) and pi(z) be polynomials defined by p*(z) := 7Z 2 + f3z + 1 and 
pi(z) := p*(%(*)-p( )p*( z ) = (1 - \j\ 2 )z + 7/3) respectively. From Theorem 4.3.2 in [8], the 
eigenvalues pi of p(z) satisfy \p\\ < 1 and | >02 1 < 1 if and only if |p(0)| < |p*(0)| and the zero of 
Pi(z), which we denote by r], satisfies < 1. The first inequality is equivalent to the inequality 

1 > I7I, and the second one \n\ < 1 is equivalent to |1 — |7|| 2 > \(3 — f3j\ 2 . 
We obtain after some works 

1 - | 7 | 2 = Ak (2 - 6k + 2k 2 - k 3 + (l - 3k - 2k 2 + 2k 3 ) cos 9 - k 3 cos 2 9) 
= 4k/(k,0), 
|1-|7|| 2 -|^-/37| 2 

= (2k sin f ) 4 (3 - 12k + 11k 2 - 2k 3 + k 4 - 2k 2 (l - k + k 2 ) cos 9 + 2k 4 cos 2 9) 
= (2Ksin|) 4 5(K,6»), 

where k = A(l — A). It is straightforward to see that /(k, 9) > and §(k, 9) > for all < k < j 
and < 9 < 2tt. Indeed, 

8 K f(K, 9) = -3(cos 9 - 1)V + 4(1 - cos 9)k - 3(2 + cos 9) < - ^ - 3 cos 9 < 0. 
Hence / is monotone decreasing with respect to k for all 9. Thus, the inequality 

f{\,9) = 39+10 cxge-cos" e >Q) 

implies that /(k, 9) > 0. Finally 

d K , K g(K, 9) = 12(2-2 cos 9 + cos 29)k 2 + 12(cos 9 - 1)k - 4 cos 9 + 22 

> o / l-cos6> a\ _ 87-96 cos 0-49 cos 20+4 cos 36> ^ n 

— u k,k9\ 2(2-2 cos 6»-cos 28) ' U > ~ 2(2-2 cos 6»-cos 20) ^ U ' 

thus, d K g(K, 9) is increasing with respect to k for all 9, and 

0^(1,0) = -lO8-12co S + c OS 20 < Q _ 

Therefore g(K, 9) is decreasing with respect to k for all 9, and we have 

g(K,9) > g{0,9) = 17 °- 26c 2 ° 5 s 6 e+cos2e > 0. 

□ 
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Figure 8: 1-D Maxwell's equations with pT = e~ = 1 (x < 0.5)and [i + = 3, e + = | (0.5 < x). 
The interface condition [H] = [E] = is imposed at the interface x = 0.5 (vertical line). The 
left plots are snap shots of the numerical solution to H at t = 0, t = 0.3, t = 0.35 and t = 0.5 
from top to bottom. The right plots are the numerical solution to E at the same time. The 
mesh size is J^. 
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